This code was executed on the results after running DEEPEst with 5000 iterations on CCP0, CCP1, CCP3, and CCP4 (the 4 studies in the current version of the paper for which DEEP estimates are important; Studies 1-4).

Load data and packages

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.3
library(bayesplot)
## Warning: package 'bayesplot' was built under R version 3.5.3
## This is bayesplot version 1.7.1
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
library(DEEPEst)
library(rstan)
## Loading required package: StanHeaders
## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For improved execution time, we recommend calling
## Sys.setenv(LOCAL_CPPFLAGS = '-march=native')
## although this causes Stan to throw an error on a few processors.
load("/Users/Kelle/Documents/DEEP Kellen/redo5000/Stan_Time_ccp1_EstHier8questions.Rdata")

c1est<-read.csv("/Users/Kelle/Documents/DEEP Kellen/redo5000/StanEstimates_Time_Parameters_ccp1_EstHier8questions.csv")

1. CCP1 (first study we ran and first in the paper)

1.1 Distribution of beta and delta

ggplot(c1est, aes(x = beta)) + geom_density()

ggplot(c1est, aes(x = delta)) + geom_density()

1.2 Average rhat across the sample

Rstan documentation says they recommend “only using the sample if R-hat is less than 1.05”. Assuming they mean average R-hat less than 1.05, we seem to pass this rule in CCP1 (barely). Antonia said rhat should be under 1, which we don’t quite pass.

#for beta
rhat(hier_time, pars = "mubeta")
## [1] 1.026951
#for delta
rhat(hier_time, pars = "mudelta_phi")
## [1] 1.011547
#for theta (the error?)
rhat(hier_time, pars = "mutheta")
## [1] 1.145379

1.3 Correlation between beta and delta

#beta and delta (annual discount factor) correlation
cor(c1est$beta,c1est$delta)
## [1] 0.2899706
#beta and daily discount rate correlation
cor(c1est$beta,c1est$r)
## [1] -0.1417037
plot(c1est$beta,c1est$delta)

1.4 Pairs plots

Antonia does this one participant at a time, so it’s not clear to me what the overall takeaway is. Seems like most pairs plots between beta and delta look okay, though beta and r looks less good (perhaps because r has big positive skew?)

First 12 participants shown below

pairs(hier_time,pars=c("beta[1]", "r[1]", "delta[1]"))

pairs(hier_time,pars=c("beta[2]", "r[2]", "delta[2]"))

pairs(hier_time,pars=c("beta[3]", "r[3]", "delta[3]"))

pairs(hier_time,pars=c("beta[4]", "r[4]", "delta[4]"))

pairs(hier_time,pars=c("beta[5]", "r[5]", "delta[5]"))

pairs(hier_time,pars=c("beta[6]", "r[6]", "delta[6]"))

pairs(hier_time,pars=c("beta[7]", "r[7]", "delta[7]"))

pairs(hier_time,pars=c("beta[8]", "r[8]", "delta[8]"))

pairs(hier_time,pars=c("beta[9]", "r[9]", "delta[9]"))

pairs(hier_time,pars=c("beta[10]", "r[10]", "delta[10]"))

pairs(hier_time,pars=c("beta[11]", "r[11]", "delta[11]"))

pairs(hier_time,pars=c("beta[12]", "r[12]", "delta[12]"))

2. CCP3 (eye-tracking and prolific, combined)

2.1 Distribution of beta and delta

load("/Users/Kelle/Documents/DEEP Kellen/redo5000/Stan_Time_ccp3_EstHier8questions.Rdata")

c3est<-read.csv("/Users/Kelle/Documents/DEEP Kellen/redo5000/StanEstimates_Time_Parameters_ccp3_EstHier8questions.csv")

ggplot(c3est, aes(x = beta)) + geom_density()

ggplot(c3est, aes(x = delta)) + geom_density()

2.2 Average rhat across the sample

Rstan documentation says they recommend “only using the sample if R-hat is less than 1.05”. Assuming they mean average R-hat less than 1.05, we seem to pass this rule in CCP3.

#for beta
rhat(hier_time, pars = "mubeta")
## [1] 1.018855
#for delta
rhat(hier_time, pars = "mudelta_phi")
## [1] 1.014837
#for theta (the error?)
rhat(hier_time, pars = "mutheta")
## [1] 1.006831

2.3 Correlation between beta and delta

#beta and delta (annual discount factor) correlation
cor(c3est$beta,c3est$delta)
## [1] 0.2901697
#beta and daily discount rate correlation
cor(c3est$beta,c3est$r)
## [1] -0.1098295
plot(c3est$beta,c3est$delta)

2.4 Pairs plots

First 12 participants shown below

pairs(hier_time,pars=c("beta[1]", "r[1]", "delta[1]"))

pairs(hier_time,pars=c("beta[2]", "r[2]", "delta[2]"))

pairs(hier_time,pars=c("beta[3]", "r[3]", "delta[3]"))

pairs(hier_time,pars=c("beta[4]", "r[4]", "delta[4]"))

pairs(hier_time,pars=c("beta[5]", "r[5]", "delta[5]"))

pairs(hier_time,pars=c("beta[6]", "r[6]", "delta[6]"))

pairs(hier_time,pars=c("beta[7]", "r[7]", "delta[7]"))

pairs(hier_time,pars=c("beta[8]", "r[8]", "delta[8]"))

pairs(hier_time,pars=c("beta[9]", "r[9]", "delta[9]"))

pairs(hier_time,pars=c("beta[10]", "r[10]", "delta[10]"))

pairs(hier_time,pars=c("beta[11]", "r[11]", "delta[11]"))

pairs(hier_time,pars=c("beta[12]", "r[12]", "delta[12]"))

3. CCP4 (mostly a replication of CCP1)

3.1 Distribution of beta and delta

load("/Users/Kelle/Documents/DEEP Kellen/redo5000/Stan_Time_ccp4_EstHier12questions.Rdata")

c4est<-read.csv("/Users/Kelle/Documents/DEEP Kellen/redo5000/StanEstimates_Time_Parameters_ccp4_EstHier12questions.csv")

ggplot(c4est, aes(x = beta)) + geom_density()

ggplot(c4est, aes(x = delta)) + geom_density()

3.2 Average rhat across the sample

Rstan documentation says they recommend “only using the sample if R-hat is less than 1.05”. Assuming they mean average R-hat less than 1.05, we seem to pass this rule in CCP3.

#for beta
rhat(hier_time, pars = "mubeta")
## [1] 1.016668
#for delta
rhat(hier_time, pars = "mudelta_phi")
## [1] 1.006234
#for theta (the error?)
rhat(hier_time, pars = "mutheta")
## [1] 1.012964

3.3 Correlation between beta and delta

#beta and delta (annual discount factor) correlation
cor(c4est$beta,c4est$delta)
## [1] 0.2832149
#beta and daily discount rate correlation
cor(c4est$beta,c4est$r)
## [1] -0.1538203
plot(c4est$beta,c4est$delta)

3.4 Pairs plots

First 12 participants shown below

pairs(hier_time,pars=c("beta[1]", "r[1]", "delta[1]"))

pairs(hier_time,pars=c("beta[2]", "r[2]", "delta[2]"))

pairs(hier_time,pars=c("beta[3]", "r[3]", "delta[3]"))

pairs(hier_time,pars=c("beta[4]", "r[4]", "delta[4]"))

pairs(hier_time,pars=c("beta[5]", "r[5]", "delta[5]"))

pairs(hier_time,pars=c("beta[6]", "r[6]", "delta[6]"))

pairs(hier_time,pars=c("beta[7]", "r[7]", "delta[7]"))

pairs(hier_time,pars=c("beta[8]", "r[8]", "delta[8]"))

pairs(hier_time,pars=c("beta[9]", "r[9]", "delta[9]"))

pairs(hier_time,pars=c("beta[10]", "r[10]", "delta[10]"))

pairs(hier_time,pars=c("beta[11]", "r[11]", "delta[11]"))

pairs(hier_time,pars=c("beta[12]", "r[12]", "delta[12]"))

4. CCP0 (used New DEEP unlike 1, 3, and 4; using cash back rewards, more ecologically valid)

4.1 Distribution of beta and delta

load("/Users/Kelle/Documents/DEEP Kellen/redo5000/Stan_Time_ccp0_EstHier16questions.Rdata")

c0est<-read.csv("/Users/Kelle/Documents/DEEP Kellen/redo5000/StanEstimates_Time_Parameters_ccp0_EstHier16questions.csv")

ggplot(c0est, aes(x = beta)) + geom_density()

ggplot(c0est, aes(x = delta)) + geom_density()

4.2 Average rhat across the sample

Rstan documentation says they recommend “only using the sample if R-hat is less than 1.05”. Assuming they mean average R-hat less than 1.05, we seem to pass this rule in CCP3.

#for beta
rhat(hier_time, pars = "mubeta")
## [1] 1.050258
#for delta
rhat(hier_time, pars = "mudelta_phi")
## [1] 1.039307
#for theta (the error?)
rhat(hier_time, pars = "mutheta")
## [1] 1.013776

4.3 Correlation between beta and delta

#beta and delta (annual discount factor) correlation
cor(c0est$beta,c0est$delta)
## [1] 0.3862682
#beta and daily discount rate correlation
cor(c0est$beta,c0est$r)
## [1] -0.21745
plot(c0est$beta,c0est$delta)

4.4 Pairs plots

First 12 participants shown below

pairs(hier_time,pars=c("beta[1]", "r[1]", "delta[1]"))

pairs(hier_time,pars=c("beta[2]", "r[2]", "delta[2]"))

pairs(hier_time,pars=c("beta[3]", "r[3]", "delta[3]"))

pairs(hier_time,pars=c("beta[4]", "r[4]", "delta[4]"))

pairs(hier_time,pars=c("beta[5]", "r[5]", "delta[5]"))

pairs(hier_time,pars=c("beta[6]", "r[6]", "delta[6]"))
## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
## Binning grid too coarse for current (small) bandwidth: consider increasing
## 'gridsize'

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
## Binning grid too coarse for current (small) bandwidth: consider increasing
## 'gridsize'

pairs(hier_time,pars=c("beta[7]", "r[7]", "delta[7]"))

pairs(hier_time,pars=c("beta[8]", "r[8]", "delta[8]"))

pairs(hier_time,pars=c("beta[9]", "r[9]", "delta[9]"))

pairs(hier_time,pars=c("beta[10]", "r[10]", "delta[10]"))

pairs(hier_time,pars=c("beta[11]", "r[11]", "delta[11]"))

pairs(hier_time,pars=c("beta[12]", "r[12]", "delta[12]"))